Modelling the dynamics of turbulent floods 
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Abstract 

Consider the dynamics of turbulent flow in rivers, estuaries and 
floods. Based on the widely used k-e model for turbulence, we use 
the techniques of centre manifold theory to derive dynamical models 
for the evolution of the water depth and of vertically averaged flow 
velocity and turbulent parameters. This new model for the shallow 
water dynamics of turbulent flow: resolves the vertical structure of 
the flow and the turbulence; includes interaction between turbulence 
and long waves; and gives a rational alternative to classical models for 
turbulent environmental flows. 
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1 Introduction 

Consider the dynamics of a turbulent flow over ground, as occurs in rivers, 
channels or floods. In such flows it is the large scale horizontal variations 
which are important; the vertical structure of velocity and turbulence may be 
expected to be determined by the local conditions of the horizontal flow. In 
this situation we may seek a model of the flow which only involves "coarse" 
depth-averaged quantities. Such models have been constructed before; for 
example, Fredsoe & Deigaard |14], pp37-39] depth-average the fc-equation 
model of turbulent flow to model the dynamics of breakers on a beach, 
whereas depth-averaged k-e equations have been used by Rastogi & Rodi 



33] to model heat and mass transfer in open channels, and by Keller & 



Rodi [23] to investigate flood plain flows. The need for such sophisticated 



models was also indicated by: Peregrine [31], p97] commenting that an em- 
pirical friction law derived from channel flow underestimates the turbulence 
in breakers and surf; and Mei p6| , p485] observing that eddy viscosities need 
to be different in and outside of the surf zone. 

However, the recent development of centre manifold theory and related 
techniques presages a much deeper understanding of the process of modelling 
nonlinear dynamics and foresees the systematic reduction of many nonlinear 
problems to an underlying low-dimensional system. For example, the process 
of depth-averaging has been shown to be deficient as a modelling paradigm 
p9| . In the context of turbulent flow, we show in §|3] how the mean motions 
may be determined by a few critical modes which have a non-trivial structure 
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in the vertical, e.g. as a first approximation the horizontal velocity u and tur- 
bulent energy density k are taken to have a cube-root dependence. Moreover, 
the amplitude of these modes and their evolution may be expressed in terms 
of depth-averaged quantities. We derive the following coupled nonlinear set 
of equations to model the turbulent, large-scale flow of water over ground, 
see (|28|): 

dv d{r)u) 

m (la) 

du vu ~ , vu 2, 

— ~ -1.030— + 0.0504 - 0.243 A W 
at r/ z r/ 2 k 

+ 0.961 # \9 - - 1.105 (lb) 

dk k^ uij^ 1 Qk 

— ~ -0.993 e- 0.0927— + (0.589 + 0.516 A)— - 1.106 u— , (lc) 
at rj 2 e rj 2 ox 

de e 2 ~ . veu 2 

— ~ -2.101— + 1.552 -3.215 AW— 
ot k kr] 2 

- 0.173 Ae^ + 0.533 A^^ - (1 + 0.735 ~\)u^- . (Id) 
ox k ox ox 

Here r\ is the water depth and it, k and e are depth averaged flow velocity, 
turbulent energy and turbulent dissipation, respectively. The other variables 
appearing are: v = C M fc 2 /e measuring the local eddy diffusivity, see (^|); and 
A = ri 2 e 2 /k 3 being proportional to the ratio of the vertical mixing time to the 
turbulent eddy time, see (^5|) . For example, in §§^T3] this model is used to 
predict the flow after a dam breaks. See in Figure |] (ppOD the formation of 
a turbulent bore rushing downstream from the dam. The turbulence in the 
bore is generally highest near the front and decays away behind as seen in 
Figure [| (p |3l"[) . This gives one example of the variations in spatial structure 
of the turbulence that underlies shallow water flows. 

Modelling turbulent flow is one of the major challenges in fluid dynamics. 
While large eddies, which can be as large as the flow domain, extract energy 
from the mean flow and feed it into turbulent motion, the eddies also act 
as a vortex stretching mechanism and transfer the energy to the smallest 
scales where viscous dissipation takes place. It is the scale at which the 
dissipation occurs which determines the rate of energy dissipation. However, 
the inflow of energy into turbulent motion is a characteristic of only the large 
scale motion. In other words, the turbulent but small scale motion is often 
dominated and determined by the large scale motion and can be treated 
as a perturbation of the mean flow. The coupling of energy transportation 
and energy dissipation with the mean flow is adequately described by widely 



Mei, Roberts Si Li, February 4, 2008 



2 The k-e model of turbulent flow 



4 



used second-order closure models. In particular, the most popular choice is 
the two equation k-e model, see for example Launder et al |[25|| , Hanjalic & 
Launder Rodi H], and Speziale [TO for reviews. 



We base our analysis upon the k-e model (§0) for the turbulence un- 
derlying the free surface of fluid in a channel, river or over a flood plain. 
The resultant model (|l|) is basically a model for the evolution of vertically 
averaged quantities; the model resolves large-scale, compared to the depth, 
dynamical structures in the horizontal. It is important to note the difference 
between models obtained by depth-averaged equations (which are known to 
be incorrect [BO] for other similar long- wave dynamics), and our model which 
is, for convenience, written in terms of depth-averaged quantities. 

In §[| we derive the "coarse" low- dimensional model ([!]) from the "fine" k- 
e equations. Despite the well recognised limitations of the k-e equations as a 
model of turbulence, we anticipate that the information retained in our coarse 
model is reasonably insensitive to deficiencies in the k-e dynamics. Further 
modelling may be done via more sophisticated Reynolds stress models for 



channel flow such as that described by Gibson & Rodi [16|. One aspect to 
note is that the model we derive has no adjustable parameters — all constants 
are determined from values established for the k-e turbulence model and its 
boundary conditions. Thus the model predictions are definitive. We describe 
some example solutions in §|5] to illustrate the dynamical predictions of the 
model. 

Finally, in Appendix [A] we comment on the status of the theory of centre 
manifolds in this development of a low-dimensional model of turbulent flow 
using centre manifold techniques. 



2 The k-e model of turbulent flow 

Consider the two-dimensional inviscid k-e model of turbulent flow over rough 
ground. Distance parallel to the ground's slope is measured by x, while we 
denote y as the direction normal to the slope. Molecular dissipation is ne- 
glected because we anticipate little direct effect for it in flood waves of a 
depth h = O(metre) over ground with roughness which may be many times 
the length-scale of viscous dissipation. Turbulent eddies are proposed to be 
the dominant mechanism for dispersion and dissipation. We denote the en- 
semble mean velocity components and pressure by u, v and p respectively, 
that is, for simplicity, we omit any distinguishing overbars (instead overbars 
will later be used to denote depth-averaged quantities). Then the incom- 
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pressible k-e model (with ensemble means) may be expressed as 





dU 

at 



du I dv_ 
dx dy 

F(p, u) 



(2a) 



where the vector u = (u, v, tj, k, e) Q is formed from the velocities u, v in the 
horizontal and vertical directions, the height of the free surface y = r)(x,t), 
the turbulent energy density k, and its dissipation rate e. The nonlinear 
model governing the evolution of the unknowns u and p is 



F{p, u) 



(2b) 
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Here the eddy viscosity 



(3) 



is a result of the turbulent mixing and varies in space and time. Later we 
use the following approximate values of the constants 



C. 



0.09 , o k = 1 



1.3, C el = 1.44, C, 



c2 



1.92. 



(4) 



in order to form definite models. Also 

2 / 



2 (tt) +2 [7r 

\dx \dy / 



f du dv " 
\dy dx 



describes the generation of turbulence through instabilities associated with 
mean-velocity gradients. T is the time-scale of the turbulent eddies and is 
frequently defined to be 



T = max jfc/e, CVy^mA 



the cut-off at viscous time-scales \Jv m /e is to avoid a singularity in turbulent 
production at a wall, see [13|, p470] for example. However, here we eschew 

1 We adopt the notation that a vector in parentheses, such as (it, v, r\, k, e), is a short- 
hand for the corresponding column vector. 
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the incorporation of direct viscous effects and so avoid this singularity by 
using T = k/e as the typical turbulent time-scale, where the overbars denote 
depth averages. The downward slope of the bottom 9 is assumed to be small 
and to have negligible variation. 

The boundary conditions on the bottom and the free surface are impor- 
tant in the details of the construction of the low- dimensional model. The 
following arguments lead to the given boundary conditions. 

• The standard condition is that, in view of the extremely low density of 
air, the pressure of the air on the fluid surface is effectively constant 
which we take to be zero without loss of generality. Thus the normal 
stress of the fluid across the free surface should vanish 



2, 2u 
p H — k 



dv du ( dv du^ 
dy dx ^ x \dx dy 



on y = r] . (5) 



3" i + 

This is only approximately true — corrections should exist of the order of 
p'r)', and similarly for other equations involving the free-surface. How- 



ever, the time-scale of gravity waves, ^£/(2ng), associated with the 
turbulent length-scale, £ oc fc 3//2 /e, should be typically much shorter 
than the turbulent eddy turn-over time, l/Vk (true for the scaling 
introduced in §§ |3.2[ ) and, as in |^T, §§2.3], we expect there to be lit- 



tle interaction between the turbulent fluctuations and the free surface 
dynamics. 

In this work we assume that the horizontal extent of the flood waves 
is small enough so that wind stress is negligible. This is in contrast to 
large-scale geophysical simulations, such as that by Arnold & Noye @ , 
where the wind stress is very important. Thus the fluid surface is free 
of tangential stress: 

/ 2 \ f du dv\ „ f dv du\ n 

( 1 -^)(fl» + sj +2 *(fl»-sJ =0 on ^"- (6) 

A wind stress could be incorporated into the model by appropriately 
replacing the zero right-hand side. 

Symmetry conditions for k and e on the free surface (see Arnold & 
Noye P| or Fredsoe & Deigaard |14|, pi 17] for examples) lead to 

dk 

— = on y = r] , (7) 
de 

and — = on y = i] . (8) 
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These assert that the energy in turbulent eddies cannot be lost or gained 
by transport through the free surface, and similarly for the turbulent 
dissipation. More sophisticated models of the free-surface effect upon 
turbulence by Gibson & Rodi Jl6|, p238] use these zero net flux bound- 
ary conditions. Spilling breakers on the water surface could perhaps 
be modelled by a turbulent production term on the right-hand side of 
these boundary conditions. 

• At the ground, y = 0, there must be no flow across the flat bottom: 

v{x,0,t) = 0. (9) 



Other boundary conditions on the ground are more arguable (com- 
pare our treatment with that of Arnold & Noye ||, ||). We are only 
interested in the flow outside of any molecular boundary layer that 
may exist on the stream bed — we imagine that the structure of any 
ordered viscous layer will be broken up by the roughness. This is sup- 



ported by recent experiments by Krogstad & Antonia |24| who show 
that roughness of a wall, even on a scale l/100th the thickness of a 
turbulent boundary layer, tends to reduce the overall anisotropy of the 
turbulence. A major limitation of the k-e model is the high level of 
anisotropy near a wall, so such a reduction in anisotropy due to rough- 
ness will favour the k-e model. Note that we treat the ground as y = 
in the mathematical model even though we imagine it to be rough. In 
effect, ensemble means are also done over all realisations of a "rough" 
bed with mean slope 9 and hence mean position y = 0. 

We suppose that the bottom inhibits the turbulence in its immediate 
neighbourhood so that the turbulent energy falls to zero: 

k = on?/ = 0. (10) 

In using the k-e model for near-wall turbulence, Durbin JTT], |12| asserts 
that dk/dy = on the wall as well. However, Figure 1 by Durbin |TT[] 
shows this latter condition is only significant for the viscous boundary 
layer — a layer we ignore due to the roughness of the ground. Instead of 
this requirement, we place the fairly weak constraint on the turbulent 
dissipation: 

de 

v— ^0 asy^O. (11) 
dy 

This is weak because v oc k 2 — > as y — > 0. In essence, this condition 
asserts that the bed does not directly act as a source or sink of turbulent 
dissipation. 
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• Although v oc k 2 — > as y — > 0, we suppose that v approaches zero 
slowly enough so that turbulence is still an effective mixing mechanism 
near the bed. Thus, the ensemble mean horizontal velocity should also 
vanish on the bed: 

u = ony = 0. (12) 

These three boundary conditions on the bed are the same as those used 
by "low- Reynolds" k-e turbulent models f3(| p64]. The difference here 
is that we do not include the near wall dependence upon local Reynolds 
numbers R t = k 2 /eu m and R y = y/ky/u m because these involve molec- 
ular viscosity v m . 

The boundary conditions are different to those we used in an earlier 
treatment of this problem |27| . The difference occurs because there we 
assumed that the stress du/dy is small near the bottom and is more 
appropriate to weakly turbulent flows. Here we seek the dynamics of 
flows with a strong level of turbulence leading to the boundary condi- 
tion ([TJ). 



3 Basis of the low-dimensional model 



In this paper we consider flows that vary "slowly" in the x and t directions. In 
this context, the meaning of "slow" is that the dynamics are slower relative 
to the vertical mixing time induced by the turbulence. In particular, the 
derivatives of the flow variables with respect to x and t are small quantities 
that can be collected with the "nonlinear" part of the equations and treated 
as perturbations. Hence we rewrite the equations as 
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C(p, u) + F{p, u, A) . 



(13) 



Treating the time and horizontal-space variations and the "nonlinear" terms 
as small, one sees that £(p, u) comprises the leading term in the equation. 

With a little adaptation, the operator £ has a critical space of equilibrium 
points parametrised by the water depth rj and the mean fields u, k and e. To 
ensure that the turbulent energy and turbulent dissipation are critical modes 
of £ and thus retained in our model of turbulent floods, we need k and e 
to be conserved to leading order. The parameter A G [0, 1] is an artificial 
parameter which we use to adjust the rate of decay of turbulent energy and 
its dissipation; by making A small we initially neglect the natural turbulent 
decay, but when A = 1 we recover the standard k-e model. This is reasonable 
because the combined effect of k = — e and e = —C t2 e 2 /k is a slow algebraic 
decay that is appropriate for the long-term centre manifold dynamics. 

3.1 Vertical mixing 

The operator £ can be considered as determining the dominant features of 
the evolution of the k-e model (0). It in turn is primarily composed of the 
differential operator 



which represents vertical mixing by turbulent eddies. By identifying this as 
the dominant term in the equations we are physically supposing that turbu- 
lent mixing is stronger than the other processes that redistribute momentum 
and turbulent energy. 

The boundary conditions on the bottom ( [LPf) and flllD in conjunction 
with the k and e components of £ in admit homogeneous solutions 
k oc y 1 / 3 and e constant. Such a cube- root profile in the vertical fits with 
arguments that the turbulence should be weaker near the bottom due to its 
constraining effects. Given such a profile, the turbulent diffusivity v oc y 2 ! 3 
and so the horizontal velocity component of £ also admits homogeneous 
solutions u oc y 1 ^ 3 . Although our long-wave model will be expressed in terms 
of depth-averaged quantities, we base the vertical structure that they measure 
to be these cube-root profiles. 

Traditionally, many theoretical approaches have assumed constant or near 
constant vertical profiles (see |b|, §§4.3.1] or J32], p670] for examples), as 
indeed we have also in an earlier treatment of this problem . However, as 
seen in experiments the horizontal velocity profile is typically curved (see f4"5j . 
Fig. 12] or 0) as is the turbulent energy (see [|14], Fig. 4.25]). A logarithmic 
profile is a well established approximation; here we work with the cube-root 
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profile as it is compatible with the k-e equations, is analytically tractable, 
and is a rough initial approximation to the logarithmic profile. 

These cube-root profiles result in downwards turbulent transport of mo- 
mentum and turbulent energy with constant flux and eventual removal from 
the fluid at the bed. In order to maintain flow at leading order we modify 
the conservative free-surface boundary conditions to supply the req- 

uisite flux at leading order, and to remove the supply at higher order. We 
replace with the boundary conditions that on y = 77 



(l-a 7 ) 



2 

Vx 



' du dv 
dy dx 



1 - a 7 ) 



dy 



du\ 


1-7 


dx J 


3r] 


drj dk 


1-7 


dx dx 


3r] 



u, (14) 
k. (15) 



where the artificial parameter 7, as for A introduced earlier, is small in the 
asymptotic scheme, but eventually will be set to 1 to recover the desired 
boundary conditions (|6|-[7|). Such manipulation of the governing equations 
has worked well in developing models of the laminar viscous flow of a thin 
fluid layer @, g2f. 



The centre manifold analysis forms a power series in 7 which needs to be 
summed for 7 = 1. The parameter a 7^ 1 is free for us to choose. Initially we 
omitted a, equivalent to choosing a = 0, but after 2 years of exploration we 
determined that an Euler transformation of the series' in 7 greatly improves 
the convergence, e.g. see Hinch |20|, Chapt. 8]. Introducing a 7^ is equivalent 
to such an Euler transformation. Another view of a 7^ is that of an over- 
relaxation parameter in an iterative scheme. In this problem it appears that 
a ~ 1/2 is a good compromise between conflicting influences and is used 
henceforth. 

From the special structure of L with the modified boundary conditions 
we deduce that there are four critical modes of interest in the long-term 
dynamics. These modes correspond to the leading-order conservation, and 
hence long-life, of fluid mass, momentum, turbulent energy and turbulent 
dissipation. These modes span the space A4 



b,v) = (g(v-y),u$(*) 



1/3 



, 0, 



H, KU V - 



E 



(16) 



where U, H, K and E are arbitrary functions of x and t. Note that the 
turbulent diffusivity, u, then also varies slowly in x and £; to leading order it 
is 

16K 2 /„.\ 2/3 
^ 9E 
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As proven in Appendix [A], the dominant operator C, linearised about 
the space of equilibrium points, has eigenvalues which are all negative (due 
to the decaying dynamics of turbulent dispersion), except for the four zero 
eigenvalues corresponding to the four conserved modes. Since all other modes 
decay exponentially quickly, the long time behaviour of the flow is determined 
by the functions U(x,t), H(x,t), K(x,t) and E(x,t). Respectively, these 
represent the vertically averaged horizontal velocity, the surface elevation, the 
vertically averaged turbulent energy, and the vertically averaged turbulent 
dissipation. In essence, we construct a "vertically averaged" model, but in u 
and k there is structure in the vertical, roughly proportional to y 1 ^ 3 , whose 
amplitude we measure by the vertical average. 



3.2 Approximating the centre manifold 

Centre manifold techniques systematically develop such a model in the ver- 
tically averaged quantities. Based on the relatively low-dimensional space of 
exponentially attractive equilibria A4o, centre manifold theory || suggests 
that the nonlinear terms "bend" Mo to a nearby manifold Ai of slow evo- 
lution. Further, Ai will similarly attract exponentially quickly all solutions 
in its vicinity; in standard formulations M. is called the centre manifold. 
Once on Ai, solutions evolve slowly according to a low-dimensional system 
of evolution equations — these evolution equations form the simplified model 
of the original dynamics. This general approach to forming low-dimensional 
models of dynamical systems is reviewed by Roberts EDI. 

In this problem, M. is parametrised by the four "amplitudes," U, H, K 
and E, which are functions of x and evolve in time. Due to the difficult 
nature of the nonlinear terms in the k-e model (0) we have to be very careful 
about these amplitudes and their derivatives. We introduce two independent 
small parameters: 5, as an amplitude scale; and i? to scale spatial derivatives. 
Then we treat 

u, k = O (8 2 ) , e = (8 3 ) , 7] = h + O (<5 2 ) , and |- = O (M) . 

(17) 

Note that in this scaling, the turbulent length-scale £ oc k 3 ! 2 je = O (1). This 
ensures that the turbulent eddies modelled by k and e are not asymptotically 
larger than the water depth; large-scale horizontal eddies are resolved by 
variations in the amplitudes of the model. Also from these scalings, the 
turbulent diffusivity v = O (k 2 /e) = O (5) and thus the time-scale of vertical 
mixing is 0(1/ 8). Consequently, we must consider horizontal scales larger 
than O (1/5), which accounts for requiring the product 8d in the scaling of 
horizontal derivatives. In standard applications of centre manifold theory, 



Mei, Roberts Si Li, February 4, 2008 



3 Basis of the low-dimensional model 



12 



we are free to scale the amplitudes in any reasonable fashion or indeed to 



treat the amplitudes as independent; as discussed in |35| a change in the 
scaling just reorders the appearance of the same set of terms in the model. 
In this application of centre manifold techniques, physical considerations and 
the non-standard nonlinearities place the constraint on the scaling that J-^ = 
O {5'd) and that 9 = O (5 Zq d). Nonetheless, we exploit usefully some of the 
capability of centre manifold techniques to treat amplitudes as independent 
variables by the above introduction of two independent small parameters, 5 
and 

Two other independent small parameters in this problem are: 7, the arti- 
ficial forcing at the free surface; and A, an artificial adjustment of turbulent 
interaction. We treat all four of these parameters as independently small. 

In terms of the field variables, we define the four amplitudes U, H, K 
and E in terms of physical quantities such that 

u = 5 2 U } (18a) 

r] = h + 5 2 H, (18b) 

k = 5 2 K, (18c) 

e = 6 3 E, (18d) 

where the overbar denotes a vertical average over the whole fluid depth at 
any x and t, for example: 

1 rv 

u = — udy . (19) 



V 

Denoting the collective amplitudes by s(x, t) = (U, H, K, E), pose the 
low-dimensional assumption that the evolution of the physical variables may 
be expressed in terms of the evolution of the four amplitudes (effectively 
equivalent to the "slaving" principle of synergetics |17| ): 



(p, u) = V(y, s) such that ^ = G(s) . (20) 

In general, we cannot find these functions V and Q exactly as this would be 
tantamount to solving exactly the original equations. Instead we determine 
asymptotic approximations in the five small parameters. 

It would be decidedly awkward to explicitly write out an asymptotic 
expansion in the five parameters. But it is also inappropriate to link their 
relative magnitudes into one parameter as we need to find relatively high- 
order in 7 but not in the others. Thus we apply an iterative algorithm 
in computer algebra to find the centre manifold and the evolution thereon 
which is based directly upon the Approximation Theorem 3 in P, Ell] and its 
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variants, as explained in detail by Roberts An outline of the procedure 
follows. 

The aim is to find the functional V and evolution Q such that the pres- 
sure, velocity and turbulence fields described by (0) form actual solutions 
of the scaled turbulent equations — this ensures fidelity between our model 
and the fluid dynamics of the k-e equations. Suppose that at some stage in 
an iterative scheme we have some approximation, V and Q . We then seek a 
correction, V' and Q', to obtain a more accurate solution to the turbulence 
equations. Substituting 

~ Os ~ 

(p, u) = V + V such that — = Q + Q' 

at 

into the scaled turbulence equations then rearranging, dropping products of 
corrections, and using a leading-order approximation wherever factors multi- 



ply corrections (see [41J for details), we obtain a system of equations for the 



corrections which is of the form 

CV' = R + Eg', (21) 

where £ = dV/ds\ s= Q is the basis for the linear subspace A4 in (|T5|), and, 
most importantly, R is the residual of the scaled k-e equations using the 
current approximation V and Q. This homo logical equation is solved by 
choosing corrections Q' to the evolution so that the right-hand side is in 
the range of C, then the correction to the fields V is determined. The 
solution is made unique by requiring that the amplitudes s have some specific 
physical meaning, here the vertical averages of the fields as in Then the 
current approximation V and Q is updated. The iteration is repeated until 
the residual of the governing equations, R, becomes zero to some order of 
error, whence the centre manifold model will be accurate to the same order 
of error (by the Approximation Theorem 3 of Carr j9]). 

A computer algebra program^ was written to perform all the necessary 
detailed algebra for this physical problem. A listing is given in Appendix |B|. 
The very important feature of this iteration scheme is that it is performed 
until the residuals of the actual governing equations are zero to some order of 
error. Thus the correctness of the results that we present here is based only 
upon the correct evaluation of the residuals and upon sufficient iterations to 
drive these to zero. The key to the correctness of the results produced by 
the computer program is the proper coding of the k-e turbulence equations. 



2 The computer algebra package REDUCE was used because of its flexible "operator" 
facility. At the time of writing, information about reduce was available from Anthony 



C. Hearn, RAND, Santa Monica, CA 90407-2138, USA. frailto : reduceQrand . org| 
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These can be seen in the computed residuals within the iterative loop of the 
program. 

4 Constructing the low-dimensional model 

As a first step in constructing a dynamical model we discard any variation 
in x and any influence of slope 9. Thus we first examine the dynamics of a 
uniform layer of turbulent fluid, k and e ^ 0, slowly decelerating, u ^ 0, due 
to turbulent drag on the bed. 

4.1 The physical fields to low-order 

As discussed in § §3.2j on the vertical mixing operator, the leading order 
approximation to the shape of the centre manifold is just the solutions to 
CV = 0. We deduce 



At higher orders in the small parameters 5, A and 7 we construct more 
refined descriptions of the fluid flow and its dynamics through the evolution of 
the amplitudes. However, we leave the influence of spatial variations through 
non-zero $ and 6 until the next section. 

By iterations of the scheme outlined in the previous section we obtain 
a basic description of the turbulence production and decay. The nonlinear 
processes and boundary condition corrections modify the cube-root profile 
and simultaneously determine the slow evolution of the amplitudes. 

It is useful to record the asymptotic expansions directly in terms of phys- 
ical quantities r/, u, k and e rather than the corresponding artificially scaled 
quantities H, U, K and E. We find the following expressions for the first 
significant modifications to the fields within the fluid, written in terms of a 
scaled vertical coordinate ( = y/t] which ranges from ( = at the bed to 
C = 1 at the fluid surface: 



<^0M)|(?) 



1/3 



1/3 




V 



0. 



(22a) 



w 



3 



u 




k 



(22b) 
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k « v (Ok+[C el aMC) + cr k v 2 (C)]u 2 + (C e2 a e + a k )v 3 (C)^i^) 

e « £ + ^1^(0^ + ^2^(04^, (22d) 

i/ « ^o(0— + [-C, x o e v x (0 + <WC)] T- 
e k 

+ [C e 2<WC)-<WC)]^- (22e) 

These expressions are correct to errors O (S 6 + A 3 + 7 3 , $) where, for exam- 
ple, a multinomial term 

^A b 7 c tf d = (V + A B + 7 C ,# D ) if 4 + TT + 7,>l, andd>L>. 
V ' J A B C ~ 

The vertical structure functions occurring in the expressions on the right- 
hand side of (22) follows. 

• For the turbulent dissipation: 

e (C) = -C 4/3 - -C 2/3 + — 

243 tA/ o 81 . 405 



512" 128" 3584' 

as shown in Figure |l[ See the effect of turbulent dissipation production, 
at a rate proportional to u 2 , through the velocity shear. Since velocity 
shear is largest near the bed, as seen in the shape of e p (C) this enhances 
turbulent dissipation e near the bed. 

However, the natural turbulent dissipation within the fluid causes a 
greater decay of turbulent dissipation near the bed, due to the smaller 
turbulent energy there, and so counters this enhancement. Being pro- 
portional to 1/0, this effect on the e-profile is greatest in weakly tur- 
bulent flows. 

For the turbulent energy density: 

M0 - l<V3 + 7 (i. c v3_i c v3 
M0 = 5< V3 -5< 2 ' 3 + ^' 3 , 

„,(<•) = -1L C V 3 + —c*' 3 - —C 1 ' 3 
vs ' 256 64 s 3584 s ' 
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Figure 1: the vertical structure of the turbulent dissipation field within the 
fluid as a function of the scaled vertical coordinate ( = y/rj: , e p (C)] 
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components of turbulent energy and velocity profile 




-0.4 -0.2 0.2 0.4 0.6 0.8 1 1.2 



Figure 2: the vertical structure functions for the turbulent energy density and 
horizontal velocity fields within the fluid as a function of the scaled vertical 

coordinate ( = y/rj: - — , v (C) (with 7 = 1);- -, 10 x t>i(C); , 

10 x v 2 (C); , 10 x v 3 (C). 



Mei, Roberts & Li, February 4, 2008 



4 Constructing the low-dimensional model 



18 



as shown in Figure |2|. Observe the cube-root structure in the vertical is 
modified to Vq((). As shown in Figure @, when 7 is set to 1 to recover 
the original boundary conditions from fll5]), the cube-root dependence 
is maintained near the bed, but is effectively flattened near the fluid 
surface to closely approximate the absence of turbulent energy flux 
through the free surface. This correction is simultaneously determined 
with a corresponding decay term in the evolution equations, as seen 
below, due to the removal of the sustaining flux. We look even closer 



at the effects of modifying the free-surface boundary conditions in §§fO. 

Also the effect of turbulent production, proportional to u 2 , through the 
velocity shear is largest near the bed; as seen in the shape of v\ (C) and 
■u 2 (C) this enhances turbulent energy k near the bed. 

However, the natural turbulent dissipation within the fluid causes a 
relatively greater decay of turbulent energy near the bed as compared 
with the body of the fluid, as seen in vz(C), and so counters this en- 
hancement. Being proportional to l/z>, this effect on the /c-profile is 
greatest in weakly turbulent flows. 

The basic cube-root structure of the horizontal velocity is modified in 
exactly the same way and for the same reasons as for the basic turbulent 
energy /c-profile. 

Modifications of the velocity profile due to the turbulent production 
and dissipation occur, but they occur primarily through the indirect 
effects of modifications to the turbulent diffusivity profile z^(C)- These 
are weak due to the subtractions in (|22b| ). 

The corresponding vertical structure of the turbulent mixing coefficient 
v is shown in Figure where the five components are 

*»(<) = fc 2 ' 3 -^, 

vs; 135 s 135 s 8505 s ' 

V2(0 = -|c + i^ + |c*, 

V JQ = iL C 2/3 _ 3 c 5/3 + 

v ^ 7 448 4 16 

V JQ = ^ C 2/3 + l C 2_3 5/3 
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components of turbulent diffusion profile 




-1.5 -1 -0.5 0.5 1 1.5 

Figure 3: the vertical structure functions for the turbulent mixing coefficient 
within the fluid as a function of the scaled vertical coordinate ( = y/rj: - 

- (right), u (Q (with 7 = 1);- 10 x ^(C); , 10 x z/ 2 (C); , 

10 x z/ 3 (C); (left), MO- 
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Simultaneously with the determination of the above fields, the solvability 
condition for the linear equations of the form fl2l|) supplies terms in the 
asymptotic, low-dimensional evolution equations for the amplitudes of the 
four critical modes. Writing these in terms of physical variables we find, 
with errors O (<5 6 + A 3 + 7 s , #): 



dt] 


~ 0, 




dt 




du 


567Z/ 


di 


rsj 

81rj 2 


dk 


~ — A 


7 


di 


.8 " 


de 
di 


9 

~ — A— C e 

8 



A eu 16 

16 1 /c 81 



3 V „ 
45^ 



uu 3 
rfk 



(23a) 
(23b) 



16(7fc 



567// - 16 



81a fc rf 



243 



7 + 



?>2C e \a e 



15cr k 



e 2 256 



where 



r -2 

243 61 rfk U ' 



A; 2 

z/(x,t) = C^ — 



-u 2 , (23c) 



(23d) 



(24) 



is a measure of the local turbulent diffusivity. These form a crude approxi- 
mation to the evolution equations for the four amplitudes of the model when 
there are no horizontal variations. 

The two fifth-order evolution equations ( |23c| - |23dj ) summarise the tur- 
bulent production and decay processes. Setting the artificial parameters 
A = 7 = 1 to approximate the dynamics of the original problem, observe: 
first, the natural decay of turbulent energy and dissipation in the body of 
the fluid; second, decay of turbulent energy with coefficient proportional to 
v/rf via turbulent mixing transporting energy to the bed; and third, the 
generation of turbulence energy and its dissipation through the shear in the 
vertical, proportional to uu 2 /rf. 

The horizontal velocity evolution (|23b| ) similarly includes terms: vu/rf 
which represents the effective drag of the bottom via turbulent mixing to 
the bed; and weak cubic, proportional to &u 3 / (r] 2 k) , and linear, proportional 
to eu/k, modification of this drag through changes in the stress tensor in 
the momentum equations (note that the coefficients are the difference of 
two terms, and that with the usual values for the parameters (|j) there is 
significant cancellation). 



The free-surface stays horizontal, equation ( 23a ), because there are no 
horizontal gradients until we look at order effects in a later subsection. 
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Table 1: terms in the series expansions in 7 of selected coefficients in the 
model for homogeneous turbulent decay. The last row is the sum of known 
terms at 7 = 1. 





in du/dt 


in dk/dt 




—uu/rj 2 


—vuX/rf 


i>u 3 /ri 2 k 


vu 2 /r\ 2 


— e 


1 





+1.03889 


+0.06542 


+0.72385 


+1.03100 


7 


+0.69136 


-0.86096 


+0.01388 


-0.13682 


-0.05665 


7 2 


+0.22783 


+0.14923 


-0.01564 


-0.02193 


+0.01602 


7 3 


+0.07479 


+0.02700 


-0.00921 


+0.00569 


+0.00250 


7 4 


+0.02448 


+0.00462 


-0.00333 


+0.00814 


+0.00027 


7 5 


+0.00801 


+0.00069 


-0.00080 


+0.00549 


-0.00002 


7 6 


+0.00262 


+0.00006 


-0.00003 


+0.00300 


-0.00003 


7 7 


-0.00086 


-0.00002 


+0.00017 


+0.00147 


-0.00001 


£ 


1.02995 


0.35950 


0.05040 


0.58892 


0.99307 



4.2 Convergence in the artificial parameters 

One limitation on the accuracy of the above model is that even within the 
k-e model of turbulence the coefficients are only approximate. This is due 
to both the modification of the free surface boundary conditions on u and k 
to (|i4 HT5f) ; and the introduction of A in (|i~3"D . Although setting 7 = 1 recovers 



the original boundary conditions (|6|-|7|) and A = 1 recovers the original k-e 
model, there is no certainty that this will give a model which is a good 
approximation to the "true" system. In essence the coefficients in the model 
are multi-variable Taylor series in 7 and in A. In this subsection we present 
evidence that these series converge for 7 = A = 1 and so we can form a 
reasonable model. 

Arbitrarily high order terms in the centre manifold expansion may be 
computed in principle. Our computer algebra program currently is limited 
by memory and time constraints to about 8th order in 7 and lower orders 
in other parameters.^] This is only attainable by the simplification of setting 
the k-e parameters to the conventional numerical values in @. By executing 
the reduce program and discarding terms O (S 6 , 7 s , A 2 , we discover more 
terms in the series in 7. Listed in Table p] are the expansions of some of the 
coefficients appearing later in the models. 

Look down the columns in the table and see that the coefficients in each 



3 In some applications |2^, ^| such routine computations can be performed to 
30th order and are used to show convincingly the convergence or otherwise of the series 
expansions. 
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sequence generally decrease by at least a factor of two. This suggests that 
the radii of convergence of the series in 7 are roughly two or more. Thus 
simply evaluating the series' at 7 = 1 is reasonably good — some are shown 
in the bottom line of the table.0 

The convergence in the parameter A is problematical because it seems to 
appear always in the combination Xri 2 e/(i>k), that is Xrj 2 e 2 /k 3 . Thus con- 
vergence depends upon the properties of the solution which is generally un- 
known beforehand. We suggest that truncating to linear terms in A forms 
an adequate approximation. It seems at least self-consistent to do this as 
later homogeneous solutions, namely (PD|-P^), show a balance for the rela- 
tively small value r) 2 e 2 /k 3 ~ 0.2. Thus the nonlinear terms in A are generally 
expected to have a negligible influence in most flows of interest. 

First, the series' are summed for 7 = 1 as discussed in the previous 
subsection. Then introducing 

~ Xr] 2 e 2 c\^ 2 /^ ^vertical mixing time 
k 3 M k/e turbulent eddy time ' 

for brevity we write the model of decay of homogeneous turbulent flow as 
follows, with errors O (<5 6 , A 2 , 

|=0, ^ (26a) 

3u ~ isu - vi$ 

— = -(1.030 + 0.359 A)— + (0.0504 -0.243 A)— (26b) 
at r\ l r\ 2 k 

= _(o.0927 + 0.993A)— + (0.589 + 0.516 A)—, (26c) 
at r/ 2 e r] 1 

— = -2.101 A— + (1.552 - 3.215 A)—-. (26d) 
at 1] 2 kr] 2 

Note that setting A = 1 to recover the original problem is just equivalent 
to using A = r] 2 e 2 /k 3 . Observe that the first terms on the right-hand sides 
of the above represent decay terms through, for example in the u and k 
equations, turbulent transport to the stream bed. The last terms in the k 
and e equations represent the production of turbulence through the velocity 
shear. 



4 Actually, the introduction of the parameter a, and the selection of a = 1/2, was 
motivated by our original series in 7 exhibiting singularities for 7 ss —1, as indicated by 
Domb-Sykes plots pp| . These singularities ruined the convergence at 7 = 1. However, an 
Euler transform of the series to accelerate convergence is precisely equivalent to choosing 
non-zero a and a few numerical experiments lead to a ~ 1/2 causing good convergence. 
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One feature of the model derived here is that it has no adjustable coef- 
ficients. All constants are derived from well known physical parameters and 
accepted constants of the k-e equations. Despite its relative complexity, the 
model has been systematically derived and the constants which appear are 
well-defined. However, there are adjustable parameters, namely the order of 
truncation of the series expansions. The model (p6|), for example, contains 
just the low-order terms in expansions in S and A. 



4.3 Dynamics of spatial structure 

The leading order effect of horizontal gradients, such as due to a sloping free 
surface, is found by computing terms of order in the asymptotic expansions. 
We describe these in this subsection. 

Dominantly, horizontal gradients affect the velocity and pressure fields. 
By computing terms to order 5 3 $ we find that the velocity fields (|22a| - [22b|) 
are modified to 

v = - ffijS* + + A 3 + 7 3 + 08) j (27a) 



dx 

u = ... + 3v 3 (C)^^ + O (S 6 + A 3 + 7 3 + $ 3 ) , (27b) 
V ox 1 



where the • • • indicate the terms on the right-hand side of ( |22 b| ) , and where 
173(C) is drawn in Figure |2|. The shape of v is required by the continuity 
equation. The modification to u asserts reasonably that at low levels of tur- 
bulence, large horizontal accelerations through decreasing depth, t] x < 0, 
cause the fluid to respond with a flatter profile through a subtraction of 173(C) 

from i>o(C) as seen m Fig ure @- 

The structure of the fields within the fluid rapidly become more compli- 
cated at higher order. We do not detail the fields any more. 

By executing the computer algebra program and discarding generated 
terms O (5 6 , 7 s , A 2 , "# 2 ), we discover first order effects of horizontal variations 
with sufficient terms in the series in 7 to sum them reliably for 7 = 1. We 
find the same production and decay terms identified in (EH) and, in addition, 
extra terms in the horizontal gradients. Using the accepted values (|4]) for the 
constants of the k-e equations, we obtain the following model with our best 
estimates of its coefficients: 

drj d(7]u) 

m ai"' (28a) 

— ~ -(1.030 + 0.359 A)— + (0.0504 - 0.243 A)-^ 
ot r] 2 r] 2 k 
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+ 



0.961 -0.019 A - (0.019 - 0.087 A) 



k 



9 



dr] 
dx 



- (1.105 + 0.104 X)u— - (0.032 - 0.056 AW — 

dx k ox 

~ u 2 de 
+ (0.025 - 0.041 A) — —, 

e Ox 



(28b) 



dk 
dt 



-0.0927 



k 3 
rfe 



VU 2 



(0.025 



0.993 e+ (0.589 + 0.516 A) , 

if 

0.011 W S -|) 



(1.106 - 0.065 A)m— - (0.030 + 0.056 X)k— 

Ox Ox 

~^uk de 
+ 0.025 - 0.060A — — , 

e Ox 



(28c) 



de 
Of 



-2.101- 



e 2 



;i.552- 3.215A) 



+ (-0.006 + 0.562 A) [0 



0.173 Ae^ 

ox 



k 



dx 



~ de 
0.735 Am-. 

<7X 



(28d) 



We expect that the coefficients in the above equations, when considered as a 
model of the k-e equations given in Section ||, are accurate as shown. Except 
for the surface equation fl28a| ), the first line in each equation is the same 

73[); subsequent lines detail the 



as in the horizontally homogeneous model 
additional terms needed to begin modelling long waves. A simpler version 
of the above model, obtained by omitting terms with small coefficients, is 
recorded in the Introduction as the model ([I]). 

Equation ( |28a| ) is an exact statement of the conservation of water and is 
not modified by any higher-order effects. To 5 3 i9, it may be written 



dr] ^du 
dt dx 







(29a) 



which is a linear description of the conservation of water. Similarly, with 
9 = and in very low levels of turbulence (v ~ 0) the horizontal momentum 
equation (|28b| ) may be written to 5 3 i9 as 



(29b) 



du dr] 
— = -0.9610—. 

dt dx 

This describes the horizontal acceleration due to slope of the fluid surface. 
These last two coupled equations form a standard description of linear wave 
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dynamics except for one remarkable feature: the effect of gravity is reduced 
by the factor 0.961. For example, this would predict that even low levels of 
turbulence reduces the phase speed of waves by about two percent. As in 
thin films of viscous fluid , the phenomenon is due to the response of the 
fluid, approximately t>o(C) shown in Figure ||, being at an angle to the forcing 
1 (either due to gravity or horizontal pressure gradients) when considered in 
the space of functions on [0,7/]. Consequently, the forcing is less effective. 
Such a depression in phase speed may be observable in the propagation of 
long- waves on turbulent flow|] 

Returning to the order 6 5 momentum equation fl28bp we note several 
interesting effects. 

• The first line contains the turbulent drag terms identified in the previ- 
ous subsection. 

• The second lines describe the effects of surface and bed slope. Within 
the square brackets: 

— the first term gives the depression of wave speed discussed above; 

— the second term very weakly enhances the phase speed correction 
in turbulent flow; 

— whereas it is difficult to ascribe one definite cause to the last 
term, coefficient modifications of the form u 2 /k are common in 
this model and reflect the relative importance of the turbulence 
on the mean flow. 

• The third and fourth lines are dominated by the nonlinear advection 
term uu x , with coefficient approximately 1.1. This coefficient is larger 
than 1 because of the shear: the maximum u(y) > u advects itself 
faster than u. This third line also shows small "cross-talk" effects in 
the advection through the u 2 (log A;) and u 2 (loge)^ terms. 

The dynamics of k and e averages are given by equations ( |28c| - |28d| ) . 

• The first lines of each equation are the same turbulent production and 
decay terms identified in the previous subsection. 

5 This modelling approach shows that where there is vertical or cross-sectional structure, 
depth-averaging or cross-sectional integration is generally unsound as a modelling tool. 
The reason is that it is the size and structure of the dynamical modes which determine 
the evolution (here approximately cube-root), and not the particular amplitudes used to 
measure the motion (here depth averages). 
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• The next line in each equation may arise from the modification of the 
turbulent production through the change of the velocity profile, seen 
in ( |27b| ), due to horizontal acceleration. 

• The remaining terms simply represent horizontal advection by the fluid 
velocity. Note that different properties are advected at different effec- 
tive speedsQ as indicated by the different coefficients of the udj dx terms: 
1.1 for u and k, and 1 for e. 

The model (|1|) reported in the introduction is a simplified version of (f28"|). 
The solutions described in the next section show that the terms neglected 
from ( P%D in writing ([]]) are relatively small, contributing at most a few 
percent in the numerical balance of the terms, and so may be neglected at 
least for initial exploration. 



5 Predictions of the new model 

In this section we investigate some of the predictions the newly derived 



model (|2~8"D might make. We look at decaying turbulence, uniform flow on 
slopes, approximation to the St Venant equations, and a dam break simula- 
tion. 



5.1 Decaying turbulence 

Homogeneous turbulence decays algebraically. If there is no slope {6 = 0), 
no variations in x, no mean flow (u = 0) and no surface waves (rj = const), 
then it is consistent to seek solutions of the model (|28|) in the form k oc t~ 2 
and e oc t~ 3 . Substituting and solving for the constants of proportionality 
the model (|28| ) predicts the turbulence decays according to 

8.97 r]H~ 2 , e~12.8r/ 2 t~ 3 , V ~ 0.565 rfr 1 , A ~ 0.227. (30) 

for large time t. The turbulence ultimately decays with the balance e ~ 
0.48 ¥/ 2 /r]. 

However, the transients towards this large time behaviour may be long. 
There are two regimes of interest characterised by large and small A compared 
to 0.227 (recall from (|25f ) that A is the ratio of the vertical mixing time to 
the turbulent eddy mixing time). 

6 Though due to the nonlinear interaction terms we should really report on the speeds 
associated with the characteristics of the equations. 
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• For small A, high turbulence k and low dissipation e, the dissipation is 
roughly constant, actually 

1 

as the turbulence decays to fl30D on a time scale of approximately 
2(r/ 2 /e) 1 / 3 . 

• For large A the vertical mixing time is relatively rapid and the turbu- 
lence decays with a different power law for some time. We find that 
e oc k 2A1 which is only a little different from fl30|). The rate of decay 
towards (|30|) is relatively slow, 

(k,e)^A(t-°- 90 ,0.90t- 1 - 90 ) , 

and forms a long lasting transient. 

The above results are for stationary water. Instead, if the water is moving 
with uniform velocity on a horizontal bed then the characteristics of the 
decaying bulk motion and turbulence are different in detail. We seek solutions 
of the model (|28| ) in the form k oc t~ 2 and e oc t~ 3 , as before, but now with 
u oc £ -1 . Substituting and solving for the constants of proportionality the 
model fl2~8|) may be rewritten as a generalised eigenvalue problem for A. It is 
then straightforward to determine the only positive solution is 

u^5.52r]t-\ k ~ 20.4 r/ 2 £" 2 , e~ 41.1 77V 3 , 

v ~ 0.91 r}H~ x , A ~ 0.199. (31) 

Numerical solutions show that there are long lasting transients of a similar 
nature to those mentioned above for stationary water. We do not elaborate 
further as this class of solutions are less likely to be of interest in applications. 

5.2 Turbulence in balance with fluid flow 

Water flowing down a slope generates turbulence that provides the drag to 
balance the gravitational forcing. Let the downward bed slope be 9 7^ 0, but 
as in the previous subsection assume are no variations in x, that is, just a 
mean flow (u ^ 0) with no surface waves (77 = const). Then it is consistent 
to seek solutions of the model (|28|) in the form u oc r/ 1 / 2 ^ 3 / 2 , k oc rj6 and 
e oc r/ 1 / 2 ^ 3 / 2 . Substituting and solving for the constants of proportionality 
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leads to a nonlinearly perturbed eigenvalue problem for A which is solved 
iteratively to give 

u « 3.16 r] 1/2 (g6) 1/2 , ifew 1.957700, e ^ 1.27 r] 1/2 (g6) 3/2 , 

z> « 0.270 ri 3/2 (g6) 1/2 , A « 0.217. (32) 

We expect this flow to be established on a time scale of the vertical mixing 
which is 

T ■ - [li 

-'mix ~ / n ' 

v V g6 

Another interesting balance occurs when we assume that the production 
of turbulence parameters, k and e, equals their dissipation through natural 
dissipation and bed drag. This leads to a reduced model in the form of 
St Venant's equation used in open channel flow. Assume that at all times 
the production and dissipation of k and e in the first lines on the right-hand 
sides of ( |28c|) and Q28d|) . That is, assume the bed slope is small enough and 
the flow is evolving slowly enough that spatial and temporal gradients are 
negligible. Then seek a balance with k oc u 2 and e oc u 3 to find 

k^ 0.224 u 2 , 0.0453 u 3 /r], v « 0.0992 rju , A « 0.184. (33) 

With this balance the momentum equation ( |28 b| ) becomes 

* = _amf + o. 94g 1.12^ _ . 018 5!^ , im 

at rj \ ox J ox r] ox 

which has exactly the same form as St Venant's equation for open channel 
flow except for the small term u 2 r) x /rj. The three significant coefficients are 
worthy of comment: the self-advection coefficient of 1.12 accounts for the 
vertical nonuniformity of the velocity profile with mean u\ the influence of 
gravity is reduced to 0.94 g because, as explained earlier for Q29bp, the re- 



sponse of the fluid flow is not constant in the vertical so some of gravitational 
forcing is not used; and lastly the bed drag u 2 /r] has coefficient 0.111 which 
is larger than typical values. However, note that such drag coefficients have 
to vary depending upon the roughness of the channel bottom as often ex- 
pressed by different roughness coefficients in Manning's law || p246,e.g.] or 
p2| , pl37,e.g.]. We surmise that the flows we describe with model (j28|) have 
strong mixing in the vertical due to strong turbulence generated by a rough 
channel bed or other extremely turbulent flows such as breaking waves or 
dam spillways. 
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5.3 Simulate a dam breaking 

One of the canonical flows of shallow water occurs after a dam breaks. Here 
we simulate such a flow and resolve the water slumping downstream and be- 
coming extremely turbulent as it does so. For simplicity we use the model ([!]) 
reported in the introduction. 

Imagine a dam at x = initially holding back water of non-dimensional 
depth Tj — 1. At time t = the dam breaks and releases the water to rush 
downstream. To avoid overly poor conditioning in the numerics we let the 
water in front of the dam be of depth rj = 0.2 (all quantities will no be 
non-dimensional so that in effect g — 1). Also to smooth the initial few time 
steps we actually set rj to a tanh profile that smoothly varied between these 
extremes such that the water slope was a maximum of 2 (rather large under 
the slowly varying assumption) at the dam. The water is assumed initially 
quiescent, u = throughout, and has a low level of turbulence, somewhat 
arbitrarily chosen to be k = 0.001. Turbulent dissipation is initially set to 
A = 0.227 so that the balance of decaying turbulence ( j30f) holds throughout. 

The model (fl|) is simply discretised on a regular grid in space-time with 
a time step of At = 1/10 and space step of Aa; = 1/6. The equations are 
discretised using a six point stencil, 3 wide in space and 2 in time, and using 
second order accurate centred differencing in both space and time. This 
leads to a sparse set of nonlinear implicit equations for the variables at the 
new time given their values at any given time. The nonlinear equations are 
solved by a few Newton's iterations requiring that the maximum norm of the 
residuals be less than 10 -4 . The model ([TJ) is purely hyperbolic, so a small 
amount of spatial diffusion is incorporated to help stabilise the simulation. 
It was physically appealing to incorporate the turbulent diffusion (uu x ) x into 
the u equation and similarly for for k and e. Such diffusion made hardly any 
difference to the simulations as a whole yet usefully avoided the generation 
of unphysical and ruinous spikes in the numerical solution. The domain of 
simulation extended from six dam heights to the left behind the dam to six 
dam heights downstream. We integrated for over a time t — 7 which is long 
enough for the disturbance to reach the ends of the computational domain 
(linear waves on the dammed water having speed 1). 

The results of this simulation are shown in Figures ^ and|^. Observe that 
when the dam breaks, the water slumps down and rushes downstream in a 
turbulent bore. The bore appears undular but may be evolving towards a se- 
ries of solitary waves — they cannot be differentiated on this time-scale. The 
turbulent structure shown in Figure [5] has interesting features. Apparently 
the energy density peaks a little behind the bore, then decays approximately 
linearly with distance. It appears that up to time t = 7 the generation of 
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(b) velocity 
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Figure 4: simulation of dam breaking showing (a) the water depth 7] and 
(b) the mean downstream velocity u. Observe the formation of a bore with 
superposed waves. 

turbulence is still significantly greater than its decay as the peak is still grow- 
ing. The turbulent dissipation and eddy viscosity behave similarly, though 
the eddy viscosity appears to peak much closer to the front of the turbu- 
lent bore. The parameter A plotted in Figure [|(d) reaffirms that relatively 
small values appear relevant to flows of interest. The peak of A at the front 
of the bore predicts that there is a lot of vertical mixing at the front, but 
less so behind the bore where A is smaller. All of the above seem physically 
reasonable. 

This model that we have derived and solved in some cases of interest 
explicitly accounts for the spatio-temporal variations of the intensity and 
broad nature of the turbulence underlying the flow of shallow water. 
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A Comments on theory in this application 

This appendix addresses the connection to the rigorous theory of centre man- 
ifolds in this application. 

It must be emphasised that throughout this paper we describe the appli- 
cation of centre manifold techniques and not the application of the rigorous 
centre manifold theory. There are two main reasons for this which we elab- 
orate on below. 

Firstly, here we construct an infinite dimensional centre manifold. At 
each point in x there are four degrees of freedom, parametrised by r/, u, k 
and e; but there are an infinitude of x positions and so there is an infinite 
number of degrees of freedom in the mathematical model. However, there is 
currently very little theory on infinite dimensional centre manifolds appear- 
ing via slowly- varying approximations fl5|, e.g.], and what there is does not 



rigorously apply here, nor does it apply to many other physically interesting 
models such as dispersion in pipes ||29|| , laminar long- wave, thin-film flows 
||39|| , and the dynamics in flow reactors ||. Instead we use the formal tech- 
niques of constructing complete low-dimensional models [lCj, [34], |35|, [36], |37| , 
techniques suggested and developed by standard applications of the theory. 
We expect that eventually theory will be developed which supports the ap- 
plication of centre manifold concepts to slowly-varying approximations. 

But there is a second obstacle to supporting this model with theory. In 
standard applications of centre manifold theory the nonlinear terms in the 
original problem are required to be smooth in the neighbourhood of the 
equilibrium under consideration (here the origin, a state of no flow and no 
turbulence). However, here many nonlinear terms are definitely not smooth; 
for example, turbulent dispersion terms such as (c>~r§^J and interaction 

terms such as C e2 y are unbounded as (u, k, e) — > 0. In rigorous applications 
of centre manifold theory, one may choose the various critical modes and 
parameters to have any set of relative orders of magnitude. The resulting 
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asymptotic expressions are the same pi], it is only the sequence in which 



the terms appear that changes with a change in relative orders of magni- 
tude. Indeed, this reflects a very desirable property of a modelling procedure, 
namely that the results are essentially independent of arbitrary human-made 
assumptions (such as order of magnitude) in the analysis. However, in this 
application the highly nonsmooth nature of the original equations means that 
in order to apply the centre manifold techniques we need to choose carefully 
the various orders of magnitudes of the variables and parameters, via fllTD- 
The aim, as in all asymptotic analyses, is to obtain a tractable and physically 
relevant leading order problem. The techniques of centre manifold theory are 
then applied, it is just that the current centre manifold theorems cannot give 
rigorous justification. 

Notwithstanding these theoretical limitations we consider the systematic 
techniques are applicable because of the attractivity of the manifold of equi- 
libria, A4q, of the linear operator C. In the spirit of centre manifold theory, 
we claim that the "small nonlinear" terms on the right-hand side of (0) just 
perturb the shape of A4 to a nearby manifold M. and perturb the evolution 
thereon. Thus our last task, and the one fulfilled in this Appendix, is to 
prove the exponential decay to .Mo- 
Linearising the k-e equations (|13|) about A4 

(po, u ) = (g(h -y),j (f ) 1/3 U{x, t), 0, H(x, t), I (f ) 1/3 K(x, t), E(x, t) 

we obtain (0, du/dt) = C(p, u) where the linear operator £ is 

|- 

ay 

d ( 8- \ ri ri (~i 8 {du 2k \ s~i d (duo k o \ 
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We seek solutions proportional to exp(At). The first thing to note is that 
we address a generalised eigen-problem 



C 



as the first row of £ comes from the continuity equation. 

Thus the first row, with the boundary condition (^), gives v = for 
any eigenvalue A. Furthermore, the third row, from the vertical momentum 
equation with the pressure boundary condition (0), then gives that p = —~k 
for any A. These considerations give no constraint on the eigenvalue A. 
Ignoring the p and v components of the eigen-problem then lead to a standard 
eigen-problem; one where the operator is in block, upper-triangular form. 
Thus consider the components in turn, starting from the last, and we show 
that all eigenvalues must be non-positive and thus A4q is attractive. 

• For any eigen-function u, if turbulent dissipation e ^ 0, then 



Ta e dy 



2/3 



dy, 



Ae 



(44) 



where 



T 



is the time-scale of cross-depth turbulent diffusion. Multiplying 
by e, Jo ■ ■ ■ dy, and integrating by parts we deduce 

A /* 4/3 /oV /3 (t) 2 ^ 
Ta e fi&dy ~ ' 

provided that the boundary conditions (|43|) (assuming e y4 oo for y — ► 
h) and 

% = ° (i 2 ^) aS V ~" °' (45) 
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holds if and only if || = 0, which 



are satisfied. The equality A 
leads to 

e = E(x, t) , 

a function independent of y. Since equation ([0]) indicates nontrivial 
solutions near y = are of the form e ~ A + By 1 / 3 , then the condi- 
tion ( f45l) is effectively equivalent to ([42]). 



Further, apply Sturm-Liouville theory to (|44]) under boundary condi- 
tions fl43|) and (f42"D . Changing the vertical variable from ?/ = hz 3 to z, 
the eigenvalue problem becomes 



dz 2 



9Tcr f Xz e 



(a form of Bessel's equation Eq.9.1.51]) with the following normal 
separate boundary conditions 

de 

— = on y = and y = h. 

dz 

Applying standard Sturm-Liouville theory, see for example Birkhoff & 
Rota 0, pp. 296] or Hartman [|T^, pp. 337ff], we see that the eigenvalues 
are discrete and must tend to infinity, 



= Ai > A 2 > 



-oo . 



Thus, linearly, solutions in the neighbourhood of the manifold A4q 
are attracted exponentially quickly to it (at a rate at least as fast as 
exp(A 2 t)). 

Sturm-Liouville theory may be also applied directly for the u and k 
components to show that any eigenvalues associated primarily with 
them are discrete. We do not record the details in the following. 

Similarly, for any eigen-function u, if turbulent dissipation e = but 
the turbulent energy k ^ 0, then 



h A ' 3 d 

Ta k dy 



2/3 



dk 



+ 



2A- 



dy 3?/ 1 / 3 



Xk. 



(46) 



Multiplying this by y 2 ^ 3 we rewrite it as 



Ta k dy 



Xy 2 ' 3 k. 
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Multiplying by k, J h • • • dy and integrating by parts we deduce 



A 



h^IoV 2 [fe (y- 1/3 k)fdy 



T(r h 



Jq y 2 / 3 k 2 dy 



< 



provided (fl"ID and 



k = o (V^ 12 ) as y 



(47) 



are satisfied. The equality A = holds here if and only if |^ [y 1 ^k^j 
0, which together with boundary condition (EDI) implies 



k = K(x,t)l(l) 



1/3 



with a function X(x, t) independent of y. Since the indicial equation 
of fl4"6| ) indicates nontrivial solutions near w = are of the form k ~ 
Ay~ 2 / 3 + By 1 / 3 , then the boundary condition ( f4~7|) is equivalent to (f4~0"D. 

The only possible eigenvalue associated with non-zero 77 is 0. 

Lastly, for any eigen-function u, iie — k — rj — but the horizontal 
velocity u^O, then 



/1 4 / 3 d 



T dy 



y 



2/3 du s 



dy, 



(48) 



and we rewrite this as 

^4/3 



-1/3 ^_ 
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1/3 



Am . 



Multiplying by u, Jq ■ ■ ■ dy and integrating by parts we deduce 



provided 



A = 
and 



h^I h y 4/3 [| (y- 1/3 u)fdy 



T 



J h u 2 dy 



< 



u = o (V /6 ) as y -> , (49) 

are satisfied. Similarly to the e and A; equations, the equality A = 
holds here if and only if |^ (jj^^uj = 0. This and the boundary 
condition (EDI) yields a unique solution 



u = U(x,t)l(l) 



1/3 
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with U(x,t) independent of y. Since the indicial equation of ( pESD indi- 
cates nontrivial solutions near y = are of the form u ~ A + By 1 / 3 , 
then the boundary condition (|^) is effectively equivalent to (|37|). 

We have proven: if the boundary conditions are satisfied, then, 

except for the four-fold eigenvalue zero whose eigenfunctions span A4o, the 
eigenvalues of C are negative and bounded away from 0. Thus we expect the 
manifold (|35D is locally attractive. Further, the time-scale of this attraction 
is the cross-depth turbulent diffusion time-scale T. 



B Computer algebra 

Here we list the reduceQ computer algebra program used to derive the long- 
wave models of turbulent flow. 



The algorithm is the iterative algorithm described in jy], adapted 
to this difficult asymptotic problem. The program refines the description 
of the centre manifold and the evolution thereon until the residual of the 
governing differential equations are driven to zero, to some asymptotic error. 
The key to the correctness of the results is then in the correct coding of the 
residuals — see inside the iterative loop. 

Note that because the thickness of the film is continuously varying in 
space and time, and because of the cube-root structure in the vertical, it is 
convenient to work with equations in terms of a scaled vertical coordinate 
z = \Jy/rj so that the free surface of the film is always z — \. However, 
the turbulence equations are not explicitly rewritten in this new coordinate 
because the computer handles all the necessary details of the transformation. 

1 COMMENT Constructs a model of turbulent 2D flow of shallow water flow 

2 on a flat slope based on the k-epsilon turbulence dynamical equations. 

3 Calculates the centre manifold & reduced dynamic system on it for the 

4 k-\epsilon model with the following boundary conditions for u, v, p, k 

5 and \epsilon: u=v=k=deps/dy=0 at y=0, dk/dn=deps/dn=0 on y=eta. Fiddle 

6 the free-surface BC to linearly force u & k at the surface. This gives 

7 roughly cube-root profile in the vertical structure of u & k, whereas 

8 eps is roughly constant. Write results in terms of z=(y/eta) " (1/3) . 

9 Here scale derivatives as ddx*del for better control. 
10 

11 Created 11/11/94, last modified 8/6/99; 
12 

13 % improve output appearance 

14 on div; off allfac; on list; on revpri; 



7 At the time of writing, information about reduce was available from Anthony C. 



Hearn, RAND, Santa Monica, CA 90407-2138, USA. nailto : reduce@rand . org 
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15 factor es ,ks , eta, us ,del ,ddx,df ,g; 
16 

17 7 maximum order of calculation in del: linear=3; non-trivial=5 

18 o:=5; 

19 7 truncate in parameters: d/dx; BC kludge; k-eps fudge 

20 let { ddx~2=>0, gamm~4=>0, lamb~2=>0>; 

21 7 let ddx=>0; gamm: =del~2*gam; lamb : =lam*del~2 ; 7 for initial results 

22 cgam:=(l-gamm*l/2) ; 7 equivalent to an Euler transform of gam series 

23 theta:=del~3*ddx*thet; 
24 

25 7 turbulence constants remove to get general formulae 

26 C_m:=9/100; C_e2 : =192/100 ; C_el : =144/100 ; s_k:=l; s_e:=13/10; 
27 

28 7 FOR ALL q SUCH THAT q>o LET del~q=0$ 

29 procedure ignore_order_gt (o) ; begin 

30 IF o=3 THEN LET del~4=0; 

31 IF o=4 THEN LET del~5=0; 

32 IF o=5 THEN LET del~6=0; 

33 IF o=6 THEN LET del~7=0; 

34 IF o=7 THEN LET del~8=0; 

35 IF o=8 THEN LET del~9=0; 

36 IF o=9 THEN LET del" 10=0; 

37 IF o=10 THEN LET del~ll=0; 

38 IF o>=ll THEN LET del~12=0; 

39 end; 
40 

41 7 amplitudes and their dependences ( eta=h+hs ) 

42 depend us,x,t; 

43 depend ks,x,t; 

44 depend es,x,t; 

45 depend hs,x,t; 

46 let {df(us,t) => gu, 

47 df (hs,t) => gh, 

48 df(ks,t) => gk, 

49 df (es,t) => ge 

50 }; 
51 

52 7 since z=(y/eta) "1/3 we need the following for d/dx, d/dt & d/dy 

53 etax:=del~2*del*ddx*df (hs,x) ; 

54 procedure df dx(a) ; 

55 begin scalar aa,bb; 

56 aa:=a*del~3*ddx; bb:=a*del*ddx; 

57 return df (bb,x) +( df(aa,eta) -z/3/eta*df (aa,z) )*df(hs,x); 

58 end; 

59 procedure dfdt(a); 

60 begin scalar aa,ugh; 

61 aa:=a*del~2*del*ddx; ugh:=if gh=0 then else gh/ (del*ddx) ; 

62 return df(a,t) +( df(aa,eta) -z/3/eta*df (aa,z) )*ugh ; 

63 end; 
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64 depend z,y; 

65 let df(z,y) => 1/ (3*eta*z~2) ; 

66 fs:={z=l}$ 
67 

68 % procedures to solve for cross-stream structures 

69 operator iav; linear iav; 
linear ise; 
linear isk; 
linear isu; 
linear isp; 
linear isv; 

=> l/(p+l) 

=> ( z"(p+2) -3/(p+5) )/(p+l)/(p+2) 
=> ( z"(p+2) -z*4/(p+5) )/(p+l)/(p+4) 
=> ( z"(p+2) -z*4/(p+5) )/(p+l)/(p+2) 
=> ( z"(p+l) -1 )/(p+l) 
=> ( z"(p+l) )/(p+l) 
/2 

z~3 -1/2 )/6 
z~3 -z*2/3 )/10 
z~3 -z*2/3 )/6 
z"2 -1 )/2 
z"2 )/2 



70 operator ise; 

71 operator isk; 

72 operator isu; 

73 operator isp; 

74 operator isv; 

75 let {iav(z~~p,z 



76 
77 
78 
79 
80 
81 
82 
83 
84 
85 
86 
87 
88 
89 
90 
91 
92 
93 



,ise(z~~p,z 
, isk(z~~p,z 
, isu(z~~p,z 
,isp(z~~p,z 
, isv(z~~p,z 
,iav(z,z) => 
,ise(z,z) => 
,isk(z,z) => 
,isu(z,z) => 
,isp(z,z) => 
,isv(z,z) => 
,iav(l,z) => 
,ise(l,z) => 
,isk(l,z) => 
,isu(l,z) => 
,isp(l,z) => 
,isv(l,z) => 

>; 

94 procedure mean(a) ; 

95 procedure solv_e(a) 

96 procedure solv_k(a) 

97 procedure solv_u(a) 

98 procedure solv_p(a) 

99 procedure solv_v(a) 



z~2 -3/5 )/2 
z~2 -z*4/5 )/4 
z~2 -z*4/5 )/2 
z -1 ) 
z ) 



3*iav(a*z~2,z) ; 

9/16*9*eta~2*s_e/anu*ise(z~2*a,z) ; 
9/16*9*eta~2*s_k/anu*isk(z~2*a,z) ; 
9/16*9*eta~2/anu*isu(z"2*a,z) ; 
3*eta*isp(z~2*a,z) ; 
3*eta*isv(z~2*a,z) ; 



100 

102 7, initial approximation for the iteration 
103 

104 cbrt:=(4/3)*z$ 

105 vu:=del~2*cbrt*us; 

106 vv:=0; 

107 vp:=g*eta*(l-z~3) +del"2*8/9*ks* (-z) ; 

108 vk:=del~2*cbrt*ks; 

109 ve:=del~3*es; 

110 vnu:=c_m*vk~2/ve; ; 

111 vmu:=c_e2*ve"2/vk; 

112 vph:=0; 
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113 

114 gu:=0; 

115 gh:=0; 

116 gk:=0; 

117 ge:=0; 
118 

11 Q V V V V V V V V V V V V V V V V W V WW v v°/w vvvvv VVVVV VVVV V VVVV V vvvv vvvvv vvvvv vvvvv vvv 
-Li> /o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o/o^ 

120 % iterate 

121 anu:=c_m*ks~2/es$ °/, scaled typical diffusion 

122 sinth : =theta-theta~3/6+theta~5/120-theta~7/5040$ 

123 costh:=l-theta~2/2+theta~4/24-theta~6/720$ 

124 repeat begin 

125 write " 

126 NEXT ITERATION 

127 "; 

128 

129 begin scalar Eqc; 

130 ignore_order_gt (o) ; 

131 7, continuity equation for v 

132 Eqc:=dfdx(vu)+df (vv,y) ; 

133 ok:=if Eqc=0 then 1 else 0; 

134 vv:=vv-solv_v(Eqc) ; 

135 end; 
136 

137 7, kinematic equation for eta 

138 gh:=SUB(fs, vv-vu*etax )/del~2; 
139 

140 7. mu equation for mu of sufficient order 

141 begin scalar Eqmu; 

142 ignore_order_gt (o+3) ; 

143 Eqmu:=(vmu*vk-C_e2*ve~2)$ 

144 ok:=if ok and(Eqmu=0) then 1 else 0; 

145 vmu : =vmu-Eqmu/ (del~2*cbrt*ks) ; 

146 end; 
147 

148 7. ph production equation of order m-1 

149 begin scalar Eqph; 

150 ignore_order_gt (o-l) ; 

151 Eqph : = (vph- (df (vu, y) +df dx ( vv) ) ~2-2*df ( vv , y) ~2-2*df dx (vu) "2) $ 

152 ok:=if ok and(Eqph=0) then 1 else 0; 

153 vph : =vph-Eqph ; 

154 end; 
155 

156 7o epsilon equation of order m+1 & BC of order m 

157 begin scalar Eqeps,BCeO,BCeh,gep; 

158 ignore_order_gt (o+l) ; 

159 Eqeps:= -dfdt(ve) -lamb*vmu +C_el* (del*es/ks) *vnu*vph 

160 +l/s_e*df (vnu*df (ve ,y) ,y) -vu*dfdx(ve)-vv*df (ve,y) 

161 +l/s_e*dfdx(vnu*dfdx(ve))$ 
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162 BCeO : =del*sub (z=0 , z~2*df (ve , y) ) $ 

163 BCeh:=del*sub(f s, df (ve,y)-etax*dfdx(ve) )$ 

164 ok:=if ok and(Eqeps=0) and(BCeO=0) and(BCeh=0) then 1 else 0; 

165 gep:=+mean(Eqeps) -160/1 17*anu/eta*BCeh; 

166 ve : =ve-solv_e ( Eqeps-gep )/del 

167 +27*eta*es/(32*ks~2*C_m*del)*BCe0*((l-z)-2-l/10) ; 

168 ge:=ge+gep/del~3; 

169 end; 
170 

171 7, nu equation of order m+1 

172 begin scalar Eqnu; 

173 ignore_order_gt (o+l) ; 

174 Eqnu:=vnu*ve-C_m*vk~2$ 

175 vnu:=vnu-Eqnu/(del~3*es) ; 

176 end; 
177 

178 7. k equation & (BC of order m-1) 

179 begin scalar Eqk , BCkO , BCkh , gkp ; 

180 ignore_order_gt (o) ; 

181 Eqk:= -dfdt(vk) -lamb*ve +vnu*vph +l/s_k*df (vnu*df (vk,y) ,y) 

182 -vu*dfdx(vk)-vv*df (vk,y) +l/s_k*df dx(vnu*df dx(vk))$ 

183 BCkO:=del*sub(z=0, vk )$ 

184 BCkh : =del* sub (f s, (df (vk,y)-etax*dfdx(vk))*cgam - (1-gamm) *vk/3/eta )$ 

185 ok:=if ok and(Eqk=0) and(BCkO=0) and(BCkh=0) then 1 else 0; 

186 gkp:=7/4*mean(z~3*Eqk) -28/9*anu/s_k/eta*BCkh; 

187 vk:=vk+solv_k( -Eqk+cbrt*gkp )/del; 

188 gk:=gk+gkp/del~2; 

189 end; 

190 7o patch up nu again because of significant changes 

191 7o nu equation of order m+1 

192 begin scalar Eqnu; 

193 ignore_order_gt (o+l) ; 

194 Eqnu:=vnu*ve-C_m*vk~2$ 

195 ok:=if ok and(Eqnu=0) then 1 else 0; 

196 vnu:=vnu-Eqnu/(del~3*es) ; 

197 end; 
198 

199 7o v equation of order m-1 for p 

200 begin scalar Eqv,BCph; 

201 ignore_order_gt (o) ; 

202 Eqv:=del*( -dfdt(vv) -df (vp,y) -g*costh -2/3*df (vk,y) 

203 -vu*dfdx(vv)-vv*df (vv,y) +dfdx( vnu*(df (vu,y)+df dx(vv) ) ) 

204 +2*df (vnu*df (vv,y) ,y) )$ 

205 BCph:=del*sub(f s, (vp+2/3*vk) * (l+etax~2) 

206 -2*vnu* (df ( vv , y ) +df dx ( vu) -etax* (df dx ( vv) +df ( vu , y) ) ) 

207 )$ 

208 ok:=if ok and(Eqv=0) and(BCph=0) then 1 else 0; 

209 vp:=vp+(solv_p(Eqv) -BCph)/del; 

210 end; 
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211 

212 % u equation for u (& BC of order m-1) 

213 begin scalar Equ , BCuO , BCuh , gup ; 

214 Equ:= -dfdt(vu) +g*sinth -dfdx(vp) -2/3*df dx(vk) 

215 -vu*df dx(vu)-vv*df (vu,y) +2*df dx(vnu*df dx(vu) ) 

216 +df (vnu*(df (vu,y)+df dx(vv) ) ,y)$ 

217 BCuO:=del*sub(z=0, vu )$ 

218 BCuh:=del*sub(f s, - (1-gamm) *vu/3/eta 

219 +cgam*((df (vu,y)+dfdx(vv))*(l-etax~2) 

220 +2*etax*(df (vv,y)-dfdx(vu))) )$ 

221 ok:=if ok and(Equ=0) and(BCuO=0) and(BCuh=0) then 1 else 0; 

222 gup:=5/4*mean(z*Equ) -20/9*anu/eta*BCuh; 

223 vu:=vu+solv_u( -Equ+gup*cbrt )/del; 

224 gu:=gu+gup/del~2; 

225 end; 
226 

227 showtime; 

228 end until ok; 
229 

230 end; 
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